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^bstracX-Tht  dipole  localization  method  (DLM)  of  the 
electroencephalogram  (EEG)  is  investigated  based  on  the  high- 
resolution  EEG.  The  finite  volume  method  (FVM)-boundary 
element  method  (BEM)  coupled  method  is  used  for  DLM.  For 
FVM,  a  novel  mesh  generation  method  is  presented  to  overcome 
the  geometric  singularity.  This  method  can  be  employed  to 
calculate  EEG  forward  problems  and  extended  to  calculate  the 
high-resolution  EEG  (HREEG).  Because  of  the  high  spatial 
resolution  of  HREEG,  as  well  as  the  merits  of  FVM  and  BEM, 
we  can  expect  that  the  FVM-BEM  coupled  method  is  more 
accurate  than  the  conventional  DLM. 
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L  INTRODUCTION 

The  investigation  of  the  brain  function  activities  is  of 
great  importance  in  the  biomedical  engineering.  The  brain 
information  with  both  high  temporal  resolution  and  high 
spatial  resolution  (SR)  is  desired.  Electroencephalogram 
(EEG)  has  an  high  temporal  resolution  within  the  millisecond 
range.  However,  the  poor  SR  of  EEG  should  be  improved. 
An  essential  way  to  enhance  SR  is  the  dipole  localization 
method  (DEM)  [1],  which  is  a  parameter  estimation 
technique  to  account  for  a  scalp  potential  distribution.  In 
recent  two  decades  the  calculation  of  the  cortical  potentials 
from  the  measured  scalp  potentials  is  developed  to  overcome 
the  blurring  effect  caused  by  the  skull,  which  is  called  high- 
resolution  EEG  (HREEG). 

For  the  realistic  models,  numerical  methods  are  necessary 
for  solving  the  EEG  problem.  The  finite  element  method 
(FEM)  and  the  boundary  element  method  (BEM)  are  used 
successfully  in  biomedical  electromagnetic  field  computation. 
BEM  is  valid  to  solve  the  models  with  homogeneous  and 
isotropic  media,  but  relatively  difficult  to  solve  the  models 
with  anisotropic  conductivity,  which  is  likely  to  be  the  reality. 
FEM  is  suitable  for  the  field  calculation,  and  it  excels  in 
handling  complicate  geometry.  However,  the  numerical 
solution  provided  by  FEM  is  not  able  to  guarantee  that  the 
current  continuity  equation  is  satisfied  precisely,  which  may 
result  in  the  calculation  error.  The  finite  volume  method 
(FVM)  is  an  alternative  approach  for  solving  EEG  problems. 
Similar  to  FEM,  FVM  is  suitable  for  handling  the 
complicated  geometry  with  anisotropic  conductivity.  The 
numerical  equations  obtained  from  the  finite  volume 
discretization  satisfy  the  current  continuity  equation.  This 
property  is  necessary  for  the  EEG  calculation,  where  the 
conductivity  may  change  greatly.  Rosenfeld  and  Tanami  [2] 
used  FVM  to  solve  the  EEG  forward  problem.  However,  they 
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adopted  hexahedron  meshes  in  which  the  geometric 
singularity  can  not  be  avoided  for  the  spheroidal  head  model. 

In  this  paper,  FVM  is  used  to  calculate  the  EEG  forward 
problem.  The  triangular  prism  meshes  are  employed  to 
overcome  the  geometric  singularity;  Then  a  layered  recursive 
algorithm  (ERA)  is  developed  to  calculate  HREEG  based  on 
FVM;  After  that,  the  FVM-BEM  coupled  approach  is 
discussed  for  DEM.  The  Eow  Resolution  Electromagnetic 
Tomography  Algorithm  (FORET A)  [3]  is  also  used  for 
reference  in  DEM.  Since  the  cortical  potential  distribution 
possesses  high  spatial  resolution  compared  with  the  scalp 
potential  distribution,  the  proposed  localization  method  based 
on  HREEG  is  expected  to  be  more  accurate  than  the 
conventional  methods. 

II.  Forward  Problem 

DEM  can  be  solved  by  an  iterative  process  of  forward 
problems.  Therefore,  the  investigation  of  efficient  methods 
for  computing  forward  problems  is  of  great  significance. 
Here  the  forward  problem  is  calculated  using  FVM. 

A.  Principle  of  FVM 

Suppose  a  dipole  source  is  positioned  in  a  conducting 
region,  the  potential  distribution  in  the  region  can  be 
expressed  by  the  current  continuity  equation,  i.e. 

=  =  /  (1) 

s  s 

where  o  is  the  conductivity  of  the  media,  /  is  the  total  current 
inside  the  closed  surfaces.  The  boundary  condition  specifies 
that  no  current  flows  into  the  air,  i.e.  adcp/dn  =  0  (^  is  the 
normal  to  the  boundary). 

To  describe  the  discretization  of  the  field  domain,  let  us 
take  a  sphere  for  instance.  The  spherical  surface  is  divided 
into  a  large  number  of  triangles.  Then  it  is  easy  to  connect  all 
the  vertices  of  the  triangles  to  the  spherical  center. 
Meanwhile  the  sphere  is  divided  into  many  concentric 
spheres.  Consequently,  a  large  number  of  triangular  prisms 
are  constructed. 

The  potential  (p  at  the  vertices  or  mesh  nodes  is  taken  as 
the  unknown  to  be  solved.  To  obtain  the  FVM  equation 
system,  a  closed  surface  around  each  node  should  be  defined 
as  the  integration  surfaces  in  (1).  First,  let  us  consider  the 
two-dimensional  case  as  shown  in  Fig.  1.  Suppose  Mi^  is  the 
point  to  be  considered,  is  the  barycentre  of  the  triangle 
Mio  Mil  Mi2.  Mioi  and  ^re  the  middle  points  of  the 
triangular  edges.  All  the  barycentres  of  triangles  around  Mi^ 
are  connected  with  the  middle  points  of  their  edges 
respectively,  so  that  a  closed  region  is  shown  hatched  in  Fig.l. 
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Now  generalize  the  two-dimensional  case  to  the  three- 
dimension  case,  as  in  Fig.  2.  Triangle  Mi^  Mi^  Mi2  is  the 
interface  of  two  adjacent  triangular  prisms.  Other  five 
triangles  surrounding  point  in  Fig.  1  can  be  obtained  by 
other  triangular  prism  interfaces.  The  hatched  region  of  Fig.  1 
is  moved  to  the  middle  section  of  the  upper  and  lower 
triangular  prisms  respectively,  which  compose  the  top  and 
the  bottom  face  of  the  closed  surface  surrounding  point 
Other  parts  of  the  closed  surfaces  are  rectangular  surfaces  as 
shown  in  Fig.  2. 

Suppose  the  number  of  the  surface  segments  of  the 
polygonal  prism  is  N.  By  taking  the  surface  of  the  prism  as 
the  integral  surface  in  (1),  we  can  get  an  equation  as 

N  N 


(2) 


here  is  the  normal  to  the  surface  segment  /.  The  integral  on 
each  surface  segment  can  be  evaluated  by  the  Gaussian 
integration  method. 

The  potential  in  a  mesh  element  is  expressed  by  the 
potentials  at  its  nodes  and  the  shape  function  of  the  triangular 
prism  as  follows 

6 

(3) 

/=1 

where  is  the  shape  function.  The  gradient  of  the  potential 
at  each  Gaussian  integration  point  on  the  surface  is  given  by 

N  dF- 
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By  substituting  (4)  to  (2),  the  discrete  equation  is  achieved, 
whose  unknowns  are  the  potentials  at  all  nodes  of  all  the 
triangular  prisms  possessing  the  considered  node.  Suppose 
the  total  number  of  nodes  is  m,  then  a  system  of  m  equations 
can  be  obtained  by  the  approach  above,  which  is  just  the 
equation  of  FVM.  By  solving  the  equation  system,  the 
potential  at  each  node  can  be  obtained. 


B.  Numerical  Simulation  of  the  Forward  Problem 

A  four-layer  sphere  model  is  used  to  represent  the  brain, 
the  cerebrospinal  fluid  (CSF),  the  skull  and  the  scalp.  The 


radius  of  the  interfaces  and  the  spherical  surface  are  6.3,  6.5, 
7.1,  and  7.5  cm  respectively.  A  unit  radial  dipole  is 
positioned  at  (5.8,  0,  0)  cm  and  its  positive  direction  is  along 
the  x-axis.  The  conductivity  of  the  cortex,  CSF  and  scalp  are 
isotropic,  which  are  0.33  S/m,  1.0  S/m,  and  0.33  S/m 
respectively.  While  the  conductivity  of  the  skull  is 
anisotropic  with  the  radial  conductivity  of  0.0042  S/m  and 
the  tangential  conductivity  of  0.042  S/m.  The  potential 
distribution  on  the  circular  boundary  formed  by  cutting  the 
sphere  along  xoz  plane  is  shown  in  Fig.  3.  The  error  of  the 
numerical  solution  compared  with  the  analytical  solution 
introduced  in  [4]  is  within  1.8%. 

III.  High-Resolution  EEG 

The  methods  for  obtaining  HREEG  are  based  on  the  fact 
that  the  region  from  the  cortex  to  the  scalp,  including  the  CSF, 
the  skull  and  the  scalp,  is  passive  electrically.  Therefore,  the 
electromagnetic  process  in  this  region  is  depicted  by  the 
Eaplace  equation.  With  the  boundary  conditions,  the  cortical 
potential  distribution  can  be  obtained. 

A.  Layered  Recursive  Algorithm  for  HREEG 

A  layered  recursive  algorithm  (ERA)  is  presented  to 
realize  HREEG.  This  algorithm  is  to  calculate  the  cortical 
potentials  from  the  known  scalp  potentials,  so  that  it  can  be 
considered  as  the  inverse  application  of  FVM. 

According  to  the  passive  property  of  the  skull,  the  right 
term  of  (2)  is  0.  Thus,  the  current  continuity  equation  for 
node  Mq  is  expressed  by 

^a>  V/  W  -0  (5) 

here  n  refers  to  the  index  of  layers  named  in  ascending  order 
from  inside  to  outside,  a^  is  the  coefficient  of  the  matrix. 

Suppose  the  node  potentials  in  layers  n+\  and  n  are  known, 
by  moving  these  terms  to  the  right  side  of  the  equation,  (5)  is 
changed  into 

(6) 

The  unknowns  in  equation  (6)  are  the  node  potentials  in  the 
layer  n-\.  (6)  is  also  expressed  as  a  matrix  equation 

According  to  the  boundary  condition  on  the  outmost  layer 


Fig.  3.  Potential  distribution  on  the  surfaee  of  a  four-layer  sphere  with 
anisotropie  eonduetivity. 
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(8) 


we  can  get  the  relationship  of  the  potential  on  the  outmost 
layer  denoted  by  n  and  that  on  the  second  layer  denoted  hy  n- 
1  as  follows 


It  is  because  that  the  side  edges  of  triangular  prisms  are  along 
the  normal  direction  of  the  surface.  Since  potential  on  layer  n 
is  known,  which  is  obtained  by  the  measurement,  the 
potential  distribution  of  the  layer  n-\  is  obtained.  From  the 
known  potential  distribution  of  layer  n  and  layer  n-l,  the 
potential  distribution  of  layer  n-2  is  obtained  by  (7). 
Consequently,  the  potential  distribution  of  each  layer  is 
deduced  recursively.  It  is  thus  called  the  “layered  recursive 
algorithm”. 


B.  Simulation  of  LRA 


A  four-layer  sphere  model  with  anisotropic  conductivity  is 
simulated.  The  radius  of  the  interfaces  and  the  spherical 
surface  are  defined  as  6.0,  6.5,  7.0,  and  7.5  cm  respectively. 
A  unit  radial  dipole  is  located  at  (5.0,  0,  0)  cm,  and  its 
positive  direction  is  along  the  x-axis.  The  conductivity  of  the 
cortex,  CSF  and  scalp  are  isotropic,  with  the  conductivity  of 
0.33  S/m,  1.0  S/m,  0.33  S/m  respectively.  The  conductivity 
of  the  skull  is  anisotropic  with  radial  conductivity  of  0.0042 
S/m  and  the  tangential  conductivity  of  0.042  S/m.  Fig.  4 
shows  the  potential  distribution  on  each  layer.  Since  the 
conductivity  in  the  skull  is  much  lower  than  the  scalp  and 
CSF,  the  potential  magnitude  decreases  rapidly.  It  means  the 
cortical  potential  is  much  greater  than  the  scalp  potential. 


IV.  Dipole  Localization 

HREEG  is  obtained  as  described  above,  which  is  of  high 
spatial  resolution  compared  with  the  potential  distribution  on 
the  scalp.  Thus,  the  dipole  localization  based  on  HREEG  can 
achieve  higher  localization  accuracy.  The  FVM-BEM 
coupled  approach  is  developed  for  the  new  DEM  based  on 
HREEG. 


Fig.  4.  Potential  distribution  on  the  layers  of  a  four-layer  sphere  with 
anisotropie  eonduetivity  by  LRA. 


A.  FVM-BEM  Coupled  Approach 


For  the  discussed  model,  FVM  is  employed  from  the  scalp 
surface  to  the  cortical  interface,  and  BEM  is  employed  inside 
the  cortical  interface. 

The  BEM  equation  is  expressed  by  [5] 

^  dV  +  Y  ff  (p— f- cos ddS  +  YB—  —  dS 
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here  (p  is  the  potential,  £  is  the  dielectric  coefficient,  r  is  a 
vector  directed  from  the  field  point  Mg  Xo  dS  ,  6  is  the  angle 
included  between  f  and  the  normal  direction  oi  dS  ,V  is  the 
field  domain,  S  is  the  boundary  of  K,  is  the  boundary 
where  the  point  Mq  is  excluded  from  V.  Since  Mq  is  on  the 
interface,  therefore  7  =  2.  According  to  the  electrostatic 
analogy  theory,  i.e.  £  is  substituded  by  a  , 
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which  is  the  potential  at  the  field  point  Y  if  the  medium  were 
homogeneous  and  of  infinite  extent,  (X,  7)  is  a  vector 

directed  from  the  source  point  X  to  the  field  point  7,  and 
dS  =  dS{X)  is  an  element  of  the  bounding  surface  S  directed 
along  the  local  normal. 

The  cortical  surface  is  divided  into  a  large  number  of  small 
triangles  the  same  as  the  discretization  of  FVM.  Hence,  the 
integral  in  (1 1)  can  be  approximated  by 

eo^OdS 


11 


<p- 


^  ff  cosGdS 

7=1 

(13) 

=e7ii 

“  dn  JJASr 

7=1 

(14) 

Fj  r  dn 

s 

here  N  is  the  total  number  of  triangles  on  the  cortical  surface. 
The  boundary  condition  on  the  cortical  interface  is 

9ci  =  9co  =  9  (15) 

=  =  (16) 
on  on 

here  and  are  the  conductivity  of  the  inside  and  outside 
of  the  interface  respectively,  (p^-  represents  the  cortical 


potential  of  the  inside  interface,  represents  the  cortical 
potential  of  the  outside  interface. 


,(«+!)  _ 


dn 


is  approximated  by 
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here  is  the  potential  on  the  cortical  interface,  is 

the  potential  just  outside  the  cortical  interface  in  the  mesh 
generation  layer,  d  is  the  distance  between  these  two  layers. 
Furthermore,  (1 1)  is  transformed  into 
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Since  the  potential  outside  the  cortical  interface  is  known, 
the  right  terms  of  the  equation  can  be  calculated,  the  only 
unknown  is  J . 


To  estimate  J,  the  Low  Resolution  Electromagnetic 
Tomography  Algorithm  (LORETA)  is  used  for  reference. 
The  principle  of  LORETA  is  equivalence  of  the  cortical 
parameter  distribution  as  the  linear  superposition  of  the 
cortical  parameter  distribution  generated  by  each  current 
dipole  source. 

Suppose  the  dipoles  in  the  brain  are  all  along  the  normal 
direction,  the  model  of  the  EEG  signal  is 

GJ  +  N  =  S  (19) 

here  S  e  mxl  is  the  right  terms  in  (18),  m  is  the  number  of 
the  vertex  on  the  cortical  interface;  J e  nxl  is  the  amplitude 
of  the  unknown  dipole  sources,  n  is  the  total  number  of  the 
discretized  points  around  the  dipole  sources;  N  is  the  noise; 
Ge  mxn  is  the  coefficient  matrix  of  the  left  term  in  (18). 
Because  it  is  only  an  analytical  calculation,  the  speed  is  high. 

The  destination  function  of  the  inverse  problem  is  the 
residual  sum  of  squares  (RSS)  between  the  measurement  of  S 
and  the  left  terms  in  (19).  It  is  solved  by  the  linear 
optimization. 

To  increase  the  stability  of  the  inverse  problem,  the 
singular  value  decomposition  (SVD)  is  used  as  following 

G  =  ULV^  (20) 

here  Z  =  diag[Xi ,  •  •  • ,  ]  is  the  singular  value,  >  0  ,  i  = 

1, ...,/;  A^- ~  0,z  =  / +  l,---,m  .  A^-  (/  =  /+!,..., m)  approximated 
to  0  is  neglected. 

In  general,  the  amplitude  of  the  dipole  on  the  true  source 
position  is  greater  than  the  amplitude  on  the  other  positions. 
Suppose  there  are  totally  n  discrete  dipole  amplitudes  to  be 
estimated,  the  dipole  with  relatively  small  amplitude  can  be 
neglected  from  the  calculation.  Therefore,  we  can  remove  a 
dipole  from  the  dipole  distribution  in  each  iterative  process. 
The  final  dipole  position  left  is  just  the  true  dipole  position. 


B.  Dipole  Localization  for  the  Realistic  Head  Model 


Fig.  5(a).  The  head  model  consists  of  four  compartments, 
which  are  the  scalp,  the  skull,  CSF  and  the  cortex  from 
outside  to  inside.  The  conductivity  are  0.33  S/m,  0.0042  S/m, 
1.0  S/m,  0.33  S/m  respectively.  One  unit  dipole  is  placed 
inside  the  cortical  interface  at  (0.00,  0.00,  6.40)  cm.  The 
potential  distribution  on  the  scalp  surface  and  the  cortical 
interface  are  shown  in  Fig.  5(b)  and  5(c)  respectively.  Based 
on  the  HREEG,  the  localization  result  is  (-0.29,  -0.10,  6.44) 
cm.  The  error  is  approximately  0.32  cm. 


(a)  (b)  (c) 

Fig.  5.  Realistic  head  model,  (a)  mesh  generation  on  the  surface  of 
the  realistic  head  model;  (b)  potential  distribution  on  the  scalp;  (c) 
potential  distribution  on  the  cortex. 


V.  Conclusions 

In  this  paper,  the  FVM-BEM  coupled  method  is  introduced 
for  DEM.  First,  a  novel  FVM  is  developed  to  solve  the 
forward  problem  and  the  HREEG.  The  geometrical 
singularity  is  overcome  by  using  the  triangular  prism  meshes. 
Then  based  on  HREEG,  the  FVM-BEM  coupled  method  is 
developed  to  localize  dipole  source.  The  simulation  results 
show  that  the  numerical  solution  is  well  accordant  with  the 
analytical  solution  in  the  forward  problem.  The  localization 
result  is  close  to  the  true  source  position  in  the  realistic  head 
model.  As  described  above,  we  can  conclude  that  FVM  is  of 
high  accuracy  in  the  realization  of  the  forward  problem  and 
HREEG.  Because  of  the  high  spatial  resolution  of  HREEG, 
as  well  as  the  merits  of  FVM  and  BEM,  the  proposed  method 
can  achieve  high  precision  in  the  dipole  localization. 
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The  approaches  above  are  extended  to  the  realistic  head 
model.  The  mesh  generation  on  the  head  surface  is  shown  in 


